Excitons in a Photosynthetic Light-Harvesting System: A Combined Molecular 
Dynamics/Quantum Chemistry and Polaron Model Study 
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The dynamics of pigment-pigment and pigment-protein interactions in light-harvesting complexes 
is studied with a novel approach that combines molecular dynamics simulations with quantum 
chemistry calculations and a polaron model analysis. The molecular dynamics simulation of light- 
harvesting complexes was performed on an 87,055 atom system comprised of an LH-II complex of 
Rhodospirillum molischianum embedded in a lipid bilayer and surrounded with appropriate water 
layers. The simulation provided information about the extent and timescales of geometrical defor- 
mations of pigment and protein residues at physiological temperatures, revealing also a pathway of 
a water molecule into the B800 binding site, as well as increased dimerization within the B850 BChl 
ring, as compared to the dimerization found for the crystal structure. For each of the 16 B850 BChls 
we performed 400 ab initio quantum chemistry calculations on geometries that emerged from the 
molecular dynamical simulations, determining the fluctuations of pigment excitation energies as a 
function of time. From the results of these calculations we construct a time-dependent Hamiltonian 
of the B850 exciton system from which we determine within linear response theory the absorp- 
tion spectrum. Finally, a polaron model is introduced to describe both the excitonic and coupled 
phonon degrees of freedom by quantum mechanics. The exciton-phonon coupling that enters into 
the polaron model, and the corresponding phonon spectral function are derived from the molecular 
dynamics and quantum chemistry simulations. The model predicts that excitons in the B850 bacte- 
riochlorophyll ring are delocalized over flve pigments at room temperature. Also, the polaron model 
permits the calculation of the absorption spectrum of the B850 excitons from the sole knowledge of 
the autocorrelation function of the excitation energies of individual BChls, which is readily available 
from the combined molecular dynamics and quantum chemistry simulations. The obtained result is 
found to be in good agreement with the experimentally measured absorption spectrum. 

PACS number(s); 87.15.Aa, 87.15. Mi, 87.16.Ac 
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I. INTRODUCTION 

Life on earth is sustained through photosynthesis. The 
initial step of photosynthesis involves absorption of light 
by the so-called light-harvesting antennae complexes, and 
funneling of the resulting electronic excitation to the pho- 
tosynthetic reaction center |^-|j] . In case of purple bac- 
teria, light-harvesting complexes are ring-shaped aggre- 
gates of proteins that contain two types of pigments, bac- 
teriochlorophylls (BChls) and carotenoids. 

Recent solution of the atomic level structure of the 
light-harvesting complex II (LH-II) of two species of pur- 
ple bacteria, Rhodopseudomonas (Rps.) acidophila Iq] 
and Rhodospirillum (Rs.) molischianum |g], has opened 
the door to an understanding of the structure- function re- 
lationship in these proteins. Fig. |l| shows the structure of 
LH-II of Rs. molischianum in its membrane-water envi- 
ronment. LH-II is an octamer, consisting of eight a- and 
eight /3- protein subunits which bind non-covalently the 
following pigments: eight BChls that absorb at 800 nm 
and are referred to as B800 BChls, 16 BChls that absorb 
at 850 nm, and are referred to as B850 BChls, and eight 
carotenoids which absorb at 500 nm. The pigment sys- 
tem of LH-II is shown in Fig. g. The difference in absorp- 
tion maxima between various pigments funnels electronic 
excitation within LH-II; pigments with higher excitation 



energy pass on their excitation to the pigments with lower 
excitation energy. As a result, energy absorbed by the 
carotenoids, or the B800 BChls, is funneled into the B850 
ring 0]. The B850 ring, in turn, transfers the electronic 
excitation to a BChl ring in another light-harvesting com- 
plex, LH-I, which is directly surrounding the photosyn- 
thetic reaction center. The B850 ring, thus, serves as a 
link in a chain of excitation transfer steps which result 
ultimately in excitation of the photosynthetic reaction 
center |^]. 

The B800 and B850 BChls have in common a tetrapy- 
rol moiety that includes the electronically excited tt- 
electron system; they differ only in the length of their 
phytyl chain B. The difference in the absorption max- 
ima arises through the difference in their ligation sites. 
Protein residues at the ligation site, especially aromatic, 
charged or polar ones, as well as hydrogen bonds be- 
tween pigment and protein residues can shift the absorp- 
tion maxima of pigments B . However, the difference in 
absorption maxima is also evoked through excitonic in- 
teractions between the B850 BChls liq]. The 7r-electron 
systems of neighboring B850 BChls are in van der Waals 
contact B , giving rise to strong couplings between their 
electronic excitations. This coupling leads to coherently 
delocalized electronic excitations, so-called excitons [ pT[ . 

What is intriguing about the light-harvesting com- 




FIG. 1. Top: Top view of the LH-II octameric complex of 
Rs. molischianum H embedded into a lipid bilayer that fits 
into a hexagonal unit cell. The simulated system consists of 
an infinite, periodic repeat of this unit cell. The lipids, 267 
molecules of palmitoyloleoylphosphatidylcholine (POPC) , are 
shown in brown, with phosphorus atoms in yellow. The a- 
and l3- protein subunits of LH-II are shown in blue and ma- 
genta, respectively; bacteriochlorophylls are shown in green, 
carotenoids in yellow. Bottom: Side view of LH-II in the 
hpid-water environment. For clarity, the front half of the 
lipids is not displayed. The color code is the same as for the 
top view; however, the a-helical segments of the protein sub- 
units are rendered as cylinders. Water molecules are shown 
in light-blue. One can recognize two ring-shaped clusters of 
BChl molecules (green). A ring of eight (only four visible) 
BChls at the top is formed by the B800 BChls; a ring of six- 
teen BChls (only partially visible) at the bottom is formed by 
16 B850 BChls. (Produced with the program VMD [S]). 
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FIG. 2. Side view of the arrangement of BChl and 
carotenoid pigments in LH-II with the protein components of 
LH-II stripped away. BChls are presented as green squares. 
Sixteen B850 BChls are arranged in the bottom ring and eight 
B800 BChls in the top ring. Eight carotenoids, shown in yel- 
low, were identified in the crystal structure at the positions 
shown (produced with the program VMD [H|). 



plexes of purple bacteria is their circular symmetry: LH- 
II of Rs. molischianum is an octamer; LH-II of Rps. aci- 
dophila is a nonamer [p|; electron density projection maps 
indicate that LH-I of Rs. rubrum is a hexadecamer [n2| . 
This poses the question whether the circular shape plays 
a functional role. In m we have constructed an effec- 
tive Hamiltonian to describe the excitons formed by the 
Qj^ excitations of B850 BChls. The Hamiltonian was as- 
sumed to exhibit a perfect 8-fold symmetry axis for a ring 
of 16 B850 BChls. In the case of such ideal symmetry, 
stationary states, the excitons, have optical properties 
characteristic of the symmetry. The dipole strength is 
unequally distributed among the excitonic states. In fact, 
two degenerate excitonic states carry almost all of the 
dipole strength. This property facilitates overall absorp- 
tion, since the two degenerate states each absorb with the 
strength of eight BChls. Furthermore, an energy trap is 
formed; the lowest exciton state carries very small dipole 
strength, making fluorescence negligible. As a result, this 
state, when populated through thermal redistribution, 
traps effectively the excitation energy. However, this sce- 
nario holds only as long as the symmetry is not broken 
significantly. 

Interactions of BChls with its protein environment in 
LH-II of -Rs. iJioUschianum, together with thermal mo- 
tions, give rise to static and dynamic disorder which man- 
ifests itself through distortions of the 8-fold symmetry. 
The classification of disorder is based on the time-scale of 
protein thermal motions compared to the characteristic 
times of absorption, fluorescence, and excitation transfer. 
If the protein motion is slow compared to the latter pro- 
cesses, the effect of the environment can be considered as 
static; otherwise, the effect is dynamic. The nature of the 
disorder (static or dynamic) is reflected in the line-shape 
functions characterizing the energy dependence of opti- 



cal absorption and emission processes. It is important to 
know these functions when studying energy transfer be- 
tween two molecules since their overlap determines the 
transfer rate. 

The disorder in LH-II affects the energy level structure 
and optical properties of excitonic states, that in turn in- 
fluence absorption, fluorescence and excitation transfer. 
To understand the dynamics of LH-II it is necessary to 
quantify the extent of disorder. One measure of the ef- 
fect of disorder is the exciton delocalization length (see 
for example |0|,|l3). For a completely symmetric aggre- 
gate at temperature T = K, the exciton delocalization 
length is equal to the total number of BChls, in the case 
of LH-II of Rs. molischianum, sixteen. With increased 
disorder the coherence of the delocalized states is lost, 
and the excitons become confined to a smaller number of 
BChls [0. 

Exciton delocalization in light-harvesting complexes 
of purple bacteria has been a subject of extensive ex- 
perimental study. Non- linear absorption pq ], circular 
dichroism M& , femtosecond transient absorption [[L7[ , flu- 
orescence p8| , pump-probe p9|] , photon echo [ po[ and 
single- molecule experiments pi|, techniques have been 
employed to study exciton delocalization length, result- 
ing in vastly different estimates which range from delo- 
calization over two BChls |g2[, to delocalization over the 
entire ring Il3] . One of the obvious explanations for dif- 
ferent estimates of exciton delocalization are differences 
in the deflnition of observables as well as differences in 
the temperature at which experiments were performed. 

It appears desirable to complement the experimental 
observations of the effect of disorder on the exciton sys- 
tem of light harvesting proteins by simulations. Such 
simulations can provide a detailed, atomic level picture 
of the disorder and its effect. However, the required simu- 
lations are extremely demanding. On the one hand, one 
needs to simulate the classical motion of light harvest- 
ing proteins in their natural lipid-water environment un- 
der periodic boundary conditions (to avoid finite size ef- 
fects) and under the infiuence of full electrostatics (dipo- 
lar forces play a key role in membrane systems). On the 
other hand, one needs to carry out, in combination with 
the classical simulations, quantum chemical calculations 
that account for the coupling between protein dynamics 
and electronic properties of the BChl systems. Fortu- 
nately, the necessary calculations have become feasible 
recently. Naturally, limitations in regard to short sam- 
pling times and still rather crude quantum chemical de- 
scriptions still exist. Nevertheless, a description of static 
and dynamic disorder based on computer simulations, 
rather than on ad hoc assumptions, would constitute a 
great step forward in our understanding of the light har- 
vesting apparatus in photosynthesis. This is particularly 
true since the apparatus is based on aggregates of nu- 
merous pigments, the functional properties of which are 
sensitive to thermal fluctuations. A model that includes 
complete information on the actual disorder and its effect 
on the electronic properties provides also an opportunity 



to improve the analysis of experimental data. 

In this paper we investigate the effect of thermal dis- 
order on excitons in LH-II at room temperature. We em- 
ploy a combined molecular dynamics/quantum chemistry 
approach to establish an effective time-dependent Hamil- 
tonian for the Qj, excitations of B850 BChls. This Hamil- 
tonian is used in a quantum mechanical treatment from 
which the spectral characteristics of the B850 BChl sys- 
tem results. The description is then recast into a polaron 
model that treats both nuclear motion and electronic de- 
grees of freedom on an equal (quantum mechanical) foot- 
ing. The theory underlying our analysis of spectral prop- 
erties is outlined in chapter II; the methods employed for 
combined molecular dynamics/quantum chemistry simu- 
lations are described in chapter HI; results, together with 
the theory of the polaron model, are presented and dis- 
cussed in chapters IV and V. Conclusions are stated in 
chapter VI. 



II. THEORY AND COMPUTATION OF 
SPECTRAL PROPERTIES 

We describe the electronic properties of the B850 
BChls in LH-II under the assumption that the rest of LH- 
II, i.e., protein, B800 BChls, and lycopenes, play merely 
the role of a thermal bath. The total Hamiltonian for 
this BChl-bath system can be written 



H ^ H 



BChl 
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BChl-Bath 



H 



Bath ■ 



(1) 



The individual Hamiltonians in Eq. {M can be further 
decomposed as 



Vr 



nuc—nuc 



HbCU — Tel + Tnuc + Vel-el + Vd-r. 
HBChl-Bath = Vel-Bath + Vnuc-Bath 
Hsath — Tsath + Vsath-Bath ■ 



The indices of the quantities introduced in Eqs.(p|]4) 
are chosen in an obvious manner to define the partial 
Hamiltonians sufficiently. Due to the large mass dif- 
ference of electrons and nuclei, we separate BChl de- 
grees of freedom into "light" and "heavy" ones, i.e., elec- 
tronic and nuclear. We introduce the following set of 
indices; r — (n.,?^2, • • • , r/v^J, for all of the electronic 
degrees of freedom of BChl; R = {Ri,R2,- ■ ■ ,Rn„^J, 
for all of the nuclear degrees of freedom of BChl; and 
Z = {Zi , Z2 , • • • , -Z^Afsatfe ) , for all of the bath degrees of 
freedom. Index Z labels bath atoms as a whole, i.e., the 
splitting of the bath degrees of freedom into nuclear and 
electronic is neglected. 

According to the Born-Oppenheimer (BO) approxima- 
tion [^ we assume that the electrons move adiabatically 
in a potential that depends parametrically on the atomic 
coordinates of BChl and bath, li Hbcm describes a single 
BChl, the above approximation applies well. The energy 
gap between ground state and lowest singlet excited state 



(2) 
(3) 
(4) 



Qj,, the relevant states, measures about 1.5 eV, which 
is too large for vibrationally induced transitions to oc- 
cur, i.e., the BO approximation holds well. However, if 
Hbcm describes an entire B850 BChl ring (comprised of 
16 BChls), energy gaps between 16 excitonic states are 
much smaller and vibrations might induce transitions be- 
tween these states. Therefore, we will assume that Hbcm 
describes only one of the B850 BChl, while the remain- 
ing B850 BChls will join the bath system. From this 
description we will derive the characteristics that enter 
into a so-called polaron model which describes vibronic 
couphng to the electronic degrees of freedom, treating 
both excitons and vibrations quantum mechanically. 

Within the framework of the BO approximation we 
construct the Schrodinger equation for the electronic 
wavefunction, which depends parametrically on the nu- 
clear coordinates R and coordinates of bath atoms Z |2^] 

HlliR,Z)^lir^;R,Z) = ej.(r,;i?,Z)C(r,;i?,Z). (5) 

Here ff*; is the electronic Hamiltonian of the i-th B850 
BChl 

HliiR,Z) = Tel + Vel-el + Vel- 



+ Vel-Bath 



(6) 

e^ and (/)^ are the a-th eigenvalue and eigenvector, re- 
spectively, of the i-th BChl. At the moment we are only 
interested in the ground and the Qy electronic state of 
BChl. The electronic problem stated in Eqs. (||,||) has 
been solved with the quantum chemistry package Gaus- 
sian 98 Ig^, as explained in section III. 

The equations of motions of the nuclear and bath de- 
grees of freedom are solved in the framework of the BO 
approximation, i.e., neglecting the non-adiabatic cou- 
pling [p3|,pi 



Hnuc+Bath{R, Z) — Tnuc + TBath + ynuc-nuc\R) (7) 
+ VBath-Bath{Z) + Vnuc-Bath{R, Z) + e"(i?, Z) . 

Dynamics of the nuclear and bath degrees of freedom 
can be described through classical molecular dynamics 
simulations. The trajectory of nuclear and bath coordi- 
nates {R{tk), Z{tk)) for various simulation times tk can 
be generated, as described also in section HI. By treat- 
ing R{tk) and Z{tk) as parameters, energies of the ground 
state, eo(tfe), and of the excited Qy state, eg {tk), will be 
calculated for each BChl i, i = 1, 2, . . . , 16. 



A. Time-dependent Hamiltonian of the B850 Ring 

We construct now the time-dependent electronic 
Hamiltonian H{t) of the BChl ring. The time- 
dependence of H{t) = H(R{t),Z{t)) arises through vi- 
brational motions of the nuclear and bath degrees of free- 
dom. The electronic Hamiltonian is expressed in the ba- 
sis comprised of the ground state of the B850 ring 



io) = ni'^o) 



(8) 



and of the sixteen Q^ excitations of individual BChls 






^i), 



(9) 



where |0q) and \4>q ) denote the ground and the Qy ex- 
cited state of the ith BChl. It is assumed that the states 
defined in (@, H) are time-independent and orthonormal, 
i.e., {i\j) = Sij for i,j — 0, 1, . . .. The Hamiltonian now 
reads [pi 



16 



H{t) = Ho{t)\0){0\ + H,{t)Y,\^{i 



where 



H,{t) - Y.^'oit) 



and 



16 

''j=i 



(10) 



(11) 



(12) 



The right-hand term in Hamiltonian dl2) describes the 
interactions between Qj, excitations of all the BChls of 
the B850 ring, and can be rewritten as 



f^iit) 



mr 



e2(i) 



\ 



W,,{t) 



W,S) 



£16 (t)/ 

(13) 



Here ei{t) = Cg [t) — €^{1) and Wi^j{t) is the induced 
dipole ~ induced dipole coupling between BChls i and j 



W,,{t) = C 



ni{t) ■ nj{t) Of i jit) ■ ni{t)) {fij{t) ■ nj{t)) 



,it)^ 



lit)' 



(14) 



where ftjit) are unit vectors describing the direction of 

the transition dipole moments dj of the ground state 
-^ Qy state transition of the j-th BChl, i.e., nj{t) = 

dj{t)/\dj{t)\, and it points from N atom of pyrol II to the 
N atom of pyrol IV in BChl j at snapshot t. rjk{t) con- 
nects the coordinates of Mg atoms of BChl j and BChl 
k at time t. C is a constant proportional to the square 
of the magnitude of the transition dipole, assumed to be 
constant, and is set to 21.12 A^ eV . The latter value 



is obtained when a computationally determined value 
of C = 64.38 A^ eV 0], which corresponds to a tran- 
sition dipole moment of H D, is rescaled to match the 
experimental value 6.3 D |^^ of the transition dipole 
moment. For simplicity we assume that the coupling 
between neighboring (and all other) BChls is given by 
Wij (t) , even though the short distance between neighbor- 
ing BChls does not justify limitation to only the leading 
term of a multipole expansion. We note, however, that 
the amplitude of calculated variations of Wij{t) is two 
orders of magnitude smaller than that for the diagonal 
matrix elements ei{t) (see results section), and that the 
error involved in the treatment will not alter the quali- 
tative behavior of the system. 

The formal solution of the time-dependent Schrodinger 
equation for Hamiltonian H{ty^'^ |23] 






at 



can be written 



u{t,to)\^ito)r' 



(15) 



(16) 



where U{t, to) is the time-development operator. By di- 
viding the time interval [t,to] into N smaller intervals, 
C/(t, io) can be written (t = tpf) 

U{t,ta) = U{tN,tN-i)UitN-i,tN-2) ■■■U{ti,ta). 

(17) 

If N is sufficiently large, H(tY^'^ can be considered to be 
constant in each time interval [i^, ifc+i], i.e., for H{tY^'^ = 
H{tkY^'' ior t G [tk,tk+i]- The time-development oper- 
ator U{tk+i, tk) can then be determined as 



U{tk+i,tk) = y^exp 



-£;W (ife+i - ife) 



(18) 



x||m)('=) W(m|| . 

Here Em and ||m)('^' are eigenvalues and eigenvectors of 
H{tkY^'^ . The eigenvectors can be expanded 



\\m) 



(fe) - 



E 



„(fc) 



(19) 



The time development operator U{t,to) can be de- 
termined from combined molecular dynamics/quantum 
chemistry data. Molecular dynamics simulation will re- 
sult in a trajectory of coordinates R{tk), Z{tk)- Based on 
the coordinates R(tk), the values of Wij{tk) in Hamil- 
tonian (13) can be calculated through Eq. (|lj). The 
diagonal matrix elements, ei(tk), are determined with 
the quantum chemistry program Gaussian 98 [E3 , as de- 
scribed in the methods section. 

The Hamiltonian H{tkY^'^ can be diagonalized, and 
the excitation energies En , and eigenvectors ||m)('') de- 
termined. Calculation of U (t,to) and |^(f))'^^'^ follows 
straightforwardly from Eqs. (p^,p7|jlq) . 



The time development operator U{t, to) for the entire 
Hamiltonian (|o|), H{t), can be easily calculated from 
U{t,to). Exploiting orthogonality between ground and 
excited states ((0|i) — 0), and noting that -ffo(i) ap- 
pearing in Eqs. (p^,pl],p^ is a c-number, one can express 



U{t,to) = exp 



J / dt'Hoit') 

" J to 

|0)(0| + Y.l^{t,t,)Uz){j\ 



(20) 



The solution of the time-dependent Schrodinger equation 
associated with H{t) is 



\^{t))=U{t,to)\^{to)). 



(21) 



B. Absorption spectrum of a coupled chlorophyll 
aggregate 

In the following we will derive an expression for the 
absorption coefficient for the B850 ring, in the framework 
of linear response theory. We note that within the dipole 
approximation the total Hamiltonian Htot describing the 
BChl system and its interaction with the radiation field 
can be chosen 



Htot = H{t) - n-Eit) 



(22) 



Here H{t) is the time-dependent Hamiltonian (HO) of the 
relevant electronic degrees of freedom of the BChl system, 
(x is the dipole moment operator that is actually a sum 
of dipole moment operators of individual BChls j, i.e.. 



E- 



(23) 



E{t) is the electric field of the monochromatic radiation 
field, conventionally given in complex form 



E{t) = Re(^oue-"^*) 



(24) 



where the unit vector u accounts for the polarization of 
the radiation field. 

The energy absorbed per unit time due to interaction 
with the field is 



dW ,'_,,, dE 



Using ( p4|) we can write 

dW i^ ^ ^ ,%,-,,/ 
— ^-Eou.{,it)){e- 



-iujt ■icji\ 



(25) 



(26) 



To calculate the absorption rate we need to determine 
the dipole moment expectation value 



^W = {u,-lit)fiUtot{t)) 



(27) into expression (pq), after some algebra, results in 



where Utot{t) is the propagator associated with the 
Hamiltonian (|22| ) governed by 

ihdtUtotit) = Htot{t)Utot{t) , (28) 

and the initial condition 

Utot{-(x^) = 1 . (29) 

We account for the effect of the radiation field in an ap- 
proximate fashion decomposing Utot{t) in the form 



Utotit) = U{t) [1 + SU{t)] 
SU{-oo) = 



(30) 
(31) 



where U{t) is defined through ( |20| ) for to = — oo. One 
can readily show that 6U{t) obeys to leading order 



ihdtSU{t) w -^(i) -Eit) 
m = U-\t)llU{t) 

The corresponding solution is 

i 



'^^(^) = ftE / dt'Mt')Ep{t') 

(3 •'-°° 



(32) 
(33) 



(34) 



where we introduced the /3-th Cartesian components of 
/2(t), E{t). From this follows, using (^, again account- 
ing only for terms to leading order 



m„(t) « Aa(i) + tJ2 dt'{[fi^{t),fip{t')])E0{t') 

n p J-oo 

(35) 

Conventionally, one writes this in the form 

/oo 
dt'x^p{t,t')Ep{t'). (36) 

where we defined the so-called susceptibility tensor 

Xa/3(t,i') = ^e(t-i')([Aa(i), A/3(i')]> (37) 

introducing the Heavyside function 0(i — t'). 

In the following we will assume that the light harvest- 
ing systems studied are initially in the electronic ground 
state |0), i.e., the susceptibility tensor is in the present 
case 



xa,/3(t,0 - -e(t-t')(o|[Aa(i),AM^')]|o). (38) 



Inserting the identity operator 



1 = |0)(0| + Y.\k){k\ 



(39) 



16 



c.0{t,t')^-Q{t-t')Y, [dkc.dip\l}{t)l}~\t')\^^ (40) 



k.l = l 



^dkpdia [U{t')U-\t) 



kl 



We have introduced here the real quantities dka = 
(0|/ia|fc) which denote the (a-th Cartesian component 
of) the transition dipole moments of the individual 
BChls k = 1,2,... 16. Using the unitarity property 
[U^Ht)]jk = [U{t)]l^ one can write (|^) 



Xat3[ 



16 



it,t')^-eit^t') J2 dkcd. 



its 



(41) 



fc,;=i 



um-'it') 



U{t)U-^it') 



kl 



from which follows 



16 



(c.f3{t,t') = --e(t-t') ^ dkcdiplm U{t)U-\t') 



k,l = l 



kl 



(42) 



Employing the notation in Eq. ( |16[ ) one can express 
lJ{t,t') = U{t)U-^{t') and conclude 



16 



Xap{t,t') 



-e(i-i') Y. dka^difjlm U{t,t') 



k,l=l 



(43) 



One can now express the dipole moment expectation 
value given by (p9) in terms of (H3). In the following we 
are actually considering the ensemble average (■ • •)e of 
the dipole moment 

/oo 
dt'{xc.f3{t,t'))eEf,{t') 
p '°° 



(44) 



,(0) 



We will assume that {rria )e is time- independent. In cal- 
culating {Xaf3{t,t'))e we are left, according to (p3[), with 

expressions of the type (d^-cd/^ U[t,t') )e. We assume 

L A kl 

that dkadifs are approximately constant over the ensem- 
ble and can be replaced by d^adip- We also assume that 
the expressions ( U{t,t') )e are translationally invari- 

L ^ kl 

ant in time, i.e., are functions of i — t' only. This permits 
us to express 

2 16 

Xc.l3{t-t') = --e(t-i') ^ dkc.di0lm{[U{t-t',O)]kl)e 



k,l = l 



k=l 



(45) 



and also 



M^{t) 



= /r^(0)\ 



E 

/3 



dt'x„p{t-t')Ep{t') 



(46) 



In the following we introduce the Fourier transform of 
the susceptibility tensor, 



Xapi^) 



dTXap{T)e'^^ 



(47) 



Since in Eq. ( p6|) both Ma{t) and £'^(i) are real, it follows 
that Yq/3(t) must also be real [cf. Eq.(42 43)], and from 
Eq. (WAj this yields the expression 



Xa/3(w) = Xl^pi-^) 



Combining (24) with ( |46| ) and ([47[), one obtains 
Af„(i) = (™i")), + i^ [xoM-'^)e^'^* 



(48) 



(49) 



We can now determine the absorption rate for radiation 



by employing this expression instead of (/t(i)) in (26). 
In doing so we consider the time average (denoted by 
an overbar) over a cycle of the radiation field and use 
exp(±iwt) = 0. This yields 
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Using Eq. (Eq) one can express 
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Expression ( pT| ) describes the absorption of radiation 
with a fixed polarization u = (1(1,^2,^3)^. We are ac- 
tually interested in absorption of light of any polariza- 
tion and coming from any direction. Carrying out the 
corresponding average (• • •)rad amounts to replacing in 
( pl| ) UaUj3 by |<5a/3. The respective absorption rate is 
described by 
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which can also be written 






(53) 



With the definition of the field intensity I — -^Eq one 
can write this 



where the absorption coefficient is [ p3| 
2Tr Lo 
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Using Eq. (Bl|) one obtains then the expression 
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(54) 
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Finally, the line shape function can be expresed in terms 
of the absorption coefficient 123] 
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where ]rf]^ is the mean square of the dipole moment of 
the system. 

The absorption coefficient for a single BChl follows 
from the above expression (p9) by setting 



Ujk{t„i,to) = Sjk JJexp 



'T^ji'^i-i) i^e ~ te-i) 



(58) 



where ei{te) is the excitation energy of BChl at time ti 
of the simulation. By inserting ( p8[ ) into Eq. (|56|), and 
using Eq. l pT\ ) one obtains a familiar expression for the 
line-shape function [Q for a single BChl 
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(59) 
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Here the sum in (^8|) has been replaced by an integral and 
eo is the average value of e{t) (we dropped the index j, 
that labels the BChls); we also defined Se{t) = e{t) — eq 
and to = 0. 



III. COMPUTATIONAL METHODS 

In this section we will describe methodological details 
of the molecular dynamics (MD) simulation on LH-II as 
well as of the quantum chemistry on BChls. 

A. Molecular Dynamics Simulation 

The simulated system, shown in Fig. ^, consists of 
the LH-II protein surrounded by a lipid bilayer and a 



35 A layer of water molecules. The system has been con- 
structed in four steps, (i) - (iv) . 

Step (i): Using X-PLOR [g8[ we added hydrogen atoms 
to the published structure of an LH-II octamer |Q , pdb 
entry ILGH. Hydrogen atoms of BChls and lycopene 
were constructed by employing the "add all hydrogens" 
command from the molecular editor in QUANTA p9[ . 
The new structure, containing explicitly all hydrogen 
atoms, was subjected to a series of minimization steps 
using X-PLOR 11 . 

Step (ii:) The modeled LH-II structure was placed into 
an elementary cell that was periodically extended. The 
cell had the shape of a hexagonal prism. Accordingly, we 
placed LH-II into a hexagonally shaped patch of lipid bi- 
layer. We chose for the hexagon a side of length a = 60 A; 
the height of the bilayer was about 42 A. 

The lipids employed were palmitoyloleoylphosphatidyl- 
choline (POPC). Ten lipid molecules were placed within 
the central opening of the LH-II (see Fig ga) and 257 
lipids were placed around the protein. The number of 
lipids in the center were chosen according to the avail- 
able area in the LH-II center and to the surface area of 
equilibrated POPC lipids [Q. The chosen number is also 
justified a posteriori since the central opening of LH-II 
remained constant in size during the subsequent equili- 
bration described below. 

Step (iii): Two water layers, with 35 A combined thick- 
ness, of the same hexagonal shape as that of the lipid bi- 
layer (a = 60 A), were added to both sides of the system, 
by removing all water molecule within 2.4 A distance 
from heavy atoms of protein or lipids. 

Step (iv): Electrostatic forces were calculated using the 
Particle-Mesh-Ewald method (PME). To assure neutral- 
ity of the system, required by this method, we added 16 
Cl~ ions to the simulated system, and placed them not 
closer than 5 A from the protein residues at positions of 
minimal electrostatic energy. However, due to the high 
mobility of the ions, the original positions chosen were 
largely irrelevant. 

Input files for the MD simulation were prepared with 
the program XPLOR ^. We used the CHARMM22 
force field M, p2| to describe protein and lipids. For 
water we employed the TIP3 model |33| , omitting in- 
ternal geometry constraints p4-pq|. Ground state par- 
tial charges for the geometry minimized structure of 
BChl (phytyl tail omitted) were calculated with the pro- 
gram JAGUAR IIj, employing the ESP method (charge, 
dipole, quadrupole and octupole moment being the ESP 
constraints) and a 6-31G** basis set. The calculated 
charges are presented in Table ||. We note here that 
we used fixed ground state charges for BChls through- 
out the simulation. The remaining force field parameters 
for the BChls were chosen as those presented in [|8[|9|. 
Force field parameters for lycopene were assigned with 
the program QUANTA ^ by exploiting QUANTA'S 
display parameter file. As pointed out above, the simu- 
lated system was first subjected to energy minimizations 
in XPLOR pq]. Additional minimizations, equilibration 
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TABLE I. Partial charges for geometry optimized 
BChl (phytyl tail omitted), calculated with the program 
JAGUAR [|||, with the ESP method (charge, dipole, 
quadrupole and octupole moment being the ESP constraints) 
and the 6-31G** basis set. To hydrogen atoms bonded to the 
same carbon we assigned an average charge. Charges for the 
phytyl tail were those presented in Pq,B9| 



and MD simulations were performed using the molecular 
dynamics program NAMD, version 2 EO]. We used the 
SHAKE constraint for all hydrogen atoms. For the di- 
electric constant we chose the value e = 1; the integration 
time-step during the equilibration phase of the simulation 
was 2 fs. The equilibration was carried out for a time pe- 
riod of 2 ns in an NpT ensemble mode with a total of 
87,055 atoms per elementary cell, at a pressure of 1 atm 
and a temperature of 300 K for periodic boundary con- 
ditions, calculating electrostatic interactions in full using 
the PME method. The cutoff for Ewald summations was 
10 A. Our choice of a 35 A water layer, and a 20 A layer of 
lipids between two neighboring LH-II molecules ensures 
adequate separation between the proteins. The choice of 
a hexagonal unit cell, as compared to a square shaped 
unit cell, significantly reduced the size of the simulated 
system. 

After the equilibration phase, we performed two molec- 
ular dynamics runs, using a 1 fs integration time-step. 
The first run was performed for 10 ps, by recording con- 
figurations every 100 fs. The second run was performed 
for 800 fs, with snapshots of the trajectory recorded ev- 
ery 2 fs. The structures associated with the latter 2 fs 
snapshots were used in quantum chemistry calculations 



of BChl's excitation energies. Due to a need to perform 
quantum chemistry calculations for all 16 BChls of the 
B850 system, and due to the high computational effort 
(calculations on four processors of Silicon Graphics Ori- 
gin 2000 required for every snapshot one hour computing 
time per BChl), we could only perform calculations for 
400 snapshots. 



B. Quantum Chemistry Calculations 

Calculations of the ground and excited state energies 
of individual BChls have been performed using the pro- 
gram Gaussian 98 [^. The coordinates and charges of 
the bath atoms provided a background charge distribu- 
tion. The calculations were performed at the Hartree- 
Fock/configuration interaction singles (HF/CIS) level, 
with the ST0-3G basis set. We used this basis set be- 
cause it is computationally the least expensive, and be- 
cause our test calculations showed that the more sophis- 
ticated 6-31G* basis set resulted in only slightly differ- 
ent fluctuations of excitation energies. To further speed 
up calculations, we restricted the active space for the 
CIS calculations to the ten highest occupied molecular 
orbitals and to the ten lowest unoccupied molecular or- 
bitals, as suggested in [ pl[ . 

The use of the ab initio (HF/CIS) method to calcu- 
late excitation energies versus semi-empirical, or classical 
calculations was discussed by the authors in [Bl|. Even 
though the ab initio method predicted incorrectly the 
average excitation energies of BChls (as opposed to the 
semi-empirical method) , it predicted well the width of the 
absorption spectrum, which is determined by the fluctua- 
tions of the energy values around the average value. We, 
thus, chose to use the ab initio method to determine the 
fluctuations, while for the average excitation energy we 
used a value of 1.57 eV. 



IV. SIMULATED PROTEIN DYNAMICS AND 
SPECTRA 

Below we will first compare the crystal structure of 
LH-II of Rs. molischianum p] , used to initiate our simu- 
lations, with the structure emerging after 2 ns MD equili- 
bration. We will then present the spectral properties pre- 
dicted for the B850 system of LH-II under the influence of 
thermal fluctuations as described in combined quantum 
chemistry, quantum mechanics, and classical mechanics 
calculations. In the following section, we will cast our 
results into a polaron model that describes the vibronic 
coupling of the B850 exciton system. 
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FIG. 3. Orientation of B850 BChls resulting from a snap- 
shot of a 10 ps long molecular dynamics simulation. Orien- 
tation of the Qy transition dipole moments is indicated by 
arrows. Angles and distances between the nearest neighbors 
(in degrees and in A), averaged over a 10 ps long simulation, 
are given on the outside and on the inside of the ring, respec- 
tively. The same numbers are shown in the right-hand side of 
the Figure, and compared to the angles and distances of the 
crystal structure. 



A. Structure of LH-II at room temperature 

In order to analyze the geometry of LH-II equilibrated 
in our MD simulations we performed a 10 ps MD simula- 
tion run. This permits us to compare the crystallograph- 
ically determined structure of LH-II of Rs. molischianum 
as reported in [p| with the simulated structure. 



1. B850 BChls 

The Mg-Mg distances between neighboring B850 
BChls, averaged over the 10 ps simulation, are shown 
in Fig. 0. The average Mg-Mg distances within and be- 
tween a/?-heterodimers are 9.8 A and 8.8 A, respectively. 
They differ from the corresponding averages in the crystal 
structure 0, which are 9.2 A within an a/3-heterodimer 
and 8.9 A between heterodimers. Average angles between 
neighboring BChls are also shown in Figure pi The aver- 
age angles between neighboring BChls are 160.9° within a 
a/3-heterodimer, and 143.5° between heterodimers. The 
corresponding angles in the crystal structure are 167.5° 
and 147.5°, respectively Q. 

We note here that the average angles and distances re- 
sult from a 10 ps sampling. The possibility that, for a 
10 ps long sampling, a structural inhomogeneity arises 
can be inferred from the fact that one of the BChls Mg- 
Mg inter-heterodimer distance was found to be much 
larger (10.4 A) than the average distance (8.6 A). This 
structural inhomogeneity will likely disappear for longer 
sampling times. 



2. B800 BChls 

The crystal structure of LH-II of Rs. molischianum 
revealed an unusual ligation of the B800 BChls [||. In 
variance with the customary ligation of BChl's central 
Mg atom to a His residue (as found for the B850 BChls), 
the Mg atom of B800 BChls is coordinated to the OSi 
atom of an Asp residue. In addition, an H2O molecule 
is found in the vicinity of the ligation site. The initial, 
minimized structure of LH-II embedded into a lipid- water 
environment did not contain any water molecules within 
the B800 binding pocket. During the course of the equili- 
bration, a water molecule diffused into the B800 binding 
pocket for seven out of eight of the B800 BChls. In Fig. 
we show a diffusion pathway of a water molecule into the 
binding site of one of the B800 BChls, recorded over the 
first 440 ps of the equilibration run. Representative dis- 
tances within the binding pocket, averaged over the 10 ps 
long simulation, are shown in Fig. 0: 2.72 A (O atom of 
H2O to OSI atom of Asp6), 3.76 A (O atom of H2O to 
Mg atom of BChl), 3.16 A (carbonyl 02 atom of BChl to 
O atom of H2O), 2.00 A (Mg atom of BChl to 0^1 atom 
of Asp6). The respective distances determined for the 
crystal structure of LH-II are 2.74 A, 4.24 A, 3.10 A and 
2.45 A, respectively M. 



B. Time-series analysis of simulation data 

The time series analysis presented in this section is 
based on an 800 fs long MD simulation. We recorded 
the nuclear (i?(ife)) and bath {Z{tk)) coordinates every 
2 fs, and generated a total of 400 snapshots. Based on 
R{tk), and Z{tk), we performed a total of 16 x 400 quan- 
tum chemistry calculations, resulting in excitation ener- 
gies ei{tk), i = 1, ..., 16, for each snapshot k. Based on 
R{tk) and Z{tk), we also calculated Wij{tk) according to 
Eq. (p^). Combining ei{tk) and Wij{tk) we constructed 
the time-dependent exciton Hamiltonian (13). 



1. Exciton Hamiltonian 

Fluctuations of two representative elements of Hamil- 
tonian matrix ([13) are presented in Fig. 0. The fluctu- 
ations of the largest off-diagonal matrix elements, i.e., 
couplings between neighboring BChls, are two orders of 
magnitude smaller than the fluctuations of diagonal ma- 
trix elements. The fluctuations of couplings between non- 
neighboring BChls are at least another order of magni- 
tude smaller than those of the nearest neighbors. This 
suggests that off-diagonal disorder is negligible compared 
to the diagonal disorder. 

Simple visual inspection of the time-dependence of e{t) 
in Fig. |5| reveals a prominent oscillatory component of 
about 20 fs. Indeed, the power spectrum of €i{t) (aver- 
aged over all 16 BChls) defined as 




FIG. 4. Binding pocket for B800 BChl and surrounding 
protein and lipid residues. Shown, in green, is a trajectory 
of a water molecule, obtained from a 440 ps long simulation. 
Snapshots of the water molecule at times t = 0(l),t = 230 ps 
(2), t = 360 ps (3), t = 362 ps (4) are represented. Distances, 
averaged over a 10 ps simulation, of the O atom of the water 
molecule, in position (4), with key protein and BChl atoms 
are given. 
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and shown in Fig. o, displays strong components in the 
range between 0.2 eV and 0.22 eV, corresponding to the 
oscillations with a period of 20.6 fs and 19.1 fs, respec- 
tively. These modes most likely originate from a C=0 or 
a methine bridge stretching vibration. We note that the 
autocorrelation function of Qj, excitation energy of BChl 
in methanol also displays a prominent vibrational mode 
of about 18 fs period pd[| . 

Fig. H reveals further vibrational modes. While the 
B850a and B850b BChls have some common modes 
(0.198 eV, 0.141 eV, 0.077 eV), some other modes peak 
at shghtly different energies; 0.218 eV, 0.21 eV, 0.182 eV, 
0.174 eV, 0.045 eV for B850a as compared to 0.214 eV, 
0.178 eV, 0.040 eV for B850b. These diflterences most 
likely stem from the difference in the protein environ- 
ment seen by B850a and B850b. 

A histogram of excitation energies efe(i) of all 16 BChls, 
is shown in Fig. |9[ Histograms of 8 B850a BChls and 8 
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FIG. 5. Matrix elements Hii{t) — ei(i) (top) and 

Hi2{t) = Wi2{t) (bottom) of Hamiltonian (|l3|), as calculated 
for the first 50 snapshots of the MD run. 
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FIG. 6. Average power spectrum of matrix elements ei(f), 
for B850a BChls i = 1, 3, . . . , 15 (full line), and for B850b 
, 16 (dashed line) . 



B850b BChls (data not shown) reveal rather large values 
for the full- width at half- maximum (FWHM), of about 
0.21 eV for B850a and 0.25 eV for B850b. Furthermore, 
the lineshape in Fig. reveals a non-Gaussian distribu- 
tion. The same non-Gaussian distribution of excitation 
energies was observed in Q|. The reason for such a dis- 
tribution might be strong coupling to the high-frequency 
modes. 

The Hamiltonian (|l^) was diagonalized for each snap- 
shot tk of the simulation, assuming that it is constant in 
the interval [tk, tfc+i]- Fluctuations of the excitonic en- 
ergies that result from the diagonalization are shown in 
Fig. 0. The Figure also shows energies of the 16 excitonic 
states averaged over the 800 fs long run, versus a corre- 
sponding average of the dipole strength. Dipolc strengths 
are given in units of individual BChl dipolc strength (1 
Sy). As can be inferred from Fig. 0, only the lowest five 
excitonic levels have dipole strengths larger than 1 . This 
suggests that in a thermally disordered exciton system 
the spectral maximum experiences a red-shift relative to 
the absorption maximum of an individual BChl. 

For comparison we show excitonic energies and dipole 
strengths for the Hamiltonian that corresponds to the 
perfect symmetric crystallographic structure. Energetic 
degeneracies of the exciton states are lifted in the fluctu- 
ating case and excitonic energies spread due to disorder. 
On the other hand, dipole strength, which is mostly con- 
tained in two degenerate excitonic levels in the symmetric 
case, is distributed more equally among the exciton states 
in case of the disordered system. The dipole forbidden 
nature of the lowest exciton state of the symmetric case 
is lifted in the disordered case. 



2. Absorption spectrum 

The absorption coefficient a{Lu) of the excitoni- 
cally coupled B850 BChls was calculated by employing 
Eqs. ( |l7| , [l8| , p6| ) . The Hamiltonian H{tk) was assumed 
to be constant in the time interval t^+i ^ tk = 0.5 fs, 
and the evolution operators U{tk+i,tk) in Eq. ( [l7| ) were 
determined accordingly. Since molecular dynamics tra- 
jectories were recorded only every 2 fs we generated ad- 
ditional data points for the Hamiltonian by linear inter- 
polation. Ensemble averaging was achieved by dividing 
the total simulation time of 800 fs into 40 fs long time 
intervals, which resulted in a total of 20 samples of 40 fs 
length. The integrand in Eq. ( pq ) is shown in Fig. ||. 
The fast fluctuation of the integrand required a choice of 
0.5 fs for the time discretization. The resulting normal- 
ized spectrum is also shown in Fig. @. 

The calculated exciton absorption spectrum is com- 
pared in Fig. @ to the experimentally measured spectrum. 
The calculated spectrum has a FWHM of about 0.138 eV, 
which is considerably larger than the FWHM of 0.05 eV 
for the experimental spectrum. Below we will discuss 
possible reasons for this discrepancy. 

We have also determined the absorption spectrum of 
individual BChls, according to Eq. (|59|). The spectrum 
was obtained from a sample of 320 (800 x 16/40) MD 
runs of 40 fs length. The integration time-step was again 
0.5 fs. The FWHM of the individual BChl absorption 
spectrum is about 0.125 eV, and its maximum is blue 
shifted with respect to the maximum of the exciton ab- 
sorption spectrum. This shift is in agreement with the 
above mentioned transfer of dipole strength into low en- 
ergy exciton states (see Fig. 0). 

Our test calculations indicate that for both the individ- 
ual BChl and exciton case, the additional peaks on both 
sides of the main absorption peak are unphysical and are 
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FIG. 7. Top, left: Excitonic energies, as determined by 
diagonalization of Hamiltonian (|l3) for the first 25 snapshots 
of the MD simulation. Energies are in eV. Bottom, left: 
Average energies of the 16 excitonic levels versus their corre- 
sponding average dipole strengths. Dipole strengths are given 
in units of dipole strength (1 Sy) of individual BChls. The 
filled circle indicates the excitation energy (1.57 eV) and the 
unit dipole strength of an individual BChl. Right: Exci- 
tonic energies for a symmetric (see text) Hamiltonian, versus 
corresponding dipole strengths. 



most likely a consequence of insufRcient sampling. 

The exciton spectrum is slightly broader than the spec- 
trum of individual BChls, apparently indicating the lack 
of an exchange narrowing of the absorption spectrum. 
Exchange narrowing is well known in J-aggregates with 
static disorder Q|, and is caused by excitonic interac- 
tions within an aggregate, i.e., the excitation experiences 
more than one realization of disorder, and, thus, effective 
disorder is an average over those realization. For purely 
statically disordered LH-II aggregates, the exchange nar- 
rowing factor was found to be about 2.8 m, in agreement 
with the theoretical value vN, where N characterizes 
the N-fold symmetry of the aggregate, i.e., in the present 
case N = 8. Exchange narrowing was also shown to occur 
for dynamically disordered aggregates. For a particular 
model, the authors in |44) have derived an expression 
for the narrowing factor which depends on the size of 
the aggregate, correlation time of the fluctuation tc, as 
well as the strength of the excitonic coupling J. For the 
values characteristic for our simulation, tc ~ 20 fs, and 
J ~ 300 cm~^, the exchange narrowing factor is about 
unity. Thus, provided that the model of Ref. Q is ap- 
plicable to our simulation, this result may explain the 
absence of exchange narrowing of our exciton absorption 
spectrum. 



3. Discussion of the results of the simulations 

For about 30 % of the MD geometries the quantum 
chemistry calculations predicted a slight degree of mixing 
between the Qj, and Q^; states. If this mixing was only an 
artifact of quantum chemistry calculations, it would have 
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FIG. 8. (a) Comparison of absorption spectra of B850 ex- 
citon (solid line) and individual BChls (dashed line), both 
obtained through Fourier transform of the time series data 
(see text) with the experimental absorption spectrum of B850 
and B800 (solid thick line) Q. (b) Time-dependence of the 
integrand in Eq. (km. 



resulted in the Qj^ energy being predicted too low, and 
fluctuation of the excitation energy to be unrealistically 
large. An immediate test of the validity of the quantum 
chemistry calculations is provided by a comparison of the 
predicted mean energy gap of 0.45 eV between the Qy 
and Qx states with the corresponding experimental value. 
Taking into account the additional shift due to excitonic 
interactions in the Qy band, of about 0.13 eV, as pre- 
dicted by our calculations, it appears that the quantum 
chemistry calculations account correctly for the experi- 
mentally measured 0.6 eV of the mean energy gap. More 
reflned quantum chemistry calculations are necessary to 
determine if the Qy — Qx mixing is indeed genuine. 

Next we discuss the choice of nearest-neighbor cou- 
plings. Various estimates of these couplings have been 
reported in the literature; for LH-II of Rs. molischianum: 
806 cm^^ and 377 cm^^ (scmicmpirical INDO/CIS cal- 
culations/effective Hamiltonian calculations [0,Q), 408 
cm~^ and 366 cm~^ (collective electronic oscillator ap- 
proach |4^); for LH-II of Rps. sphaeroides: 300 cm~^ 
and 233 cm~^ (circular dichroism studies |Q); for LH- 
II of Rps. acidophila: 238 cm"! and 213 cm^^ p!], 320 



cm and 255 cm m 
394 cm~^ and 317 cm" 



I (CIS Gaussian 94 calculations), 
(QCFF/PI calculation M), 622 
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cm-1 and 562 cm-i (ZINDO/S method ^). The aver- 
age values of couphngs used in this paper are 364 cm~^ 
and 305 cm~^. If larger values were used, as might have 
been suggested from some of the above mentioned esti- 
mates, a narrower absorption spectrum would have been 
predicted, due to an increase in the exchange narrowing 
factor. 

To reduce computational cost of our simulation, we em- 
ployed the dipole-dipole approximation to describe vari- 
ations of the nearest neighbor couplings as a function 
of time. It is well known that this approximation does 
not hold in the case when center-to-center distances be- 
tween pigments are smaller than the size of the pigments 
themselves. However, we find that the fluctuations of 
the amplitude of off-diagonal matrix elements is negligi- 
ble compared to the fluctuations of the diagonal matrix 
elements, and we expect that this would still hold if the 
proper coupling between the neighboring pigments was 
determined. 



V. ANALYSIS IN THE FRAMEWORK OF THE 
POLARON MODEL 

In order to describe the influence of dynamical disor- 
der (thermal fluctuations) on the electronic excitations of 
the B850 BChls, we will employ now a model that focuses 
on the effect of the strongest fluctuations in the exciton 
Hamiltonian (ttSh, the variation of the site energies ej{t). 
These variations arc foremost due to the coupling of the 
high-frequency vibrations to Qy excitations in BChl (see 
Fig. |6|). Initially, we will model this coupling through 
a single harmonic oscillator mode wq, that is coupled to 
each of the local BChl excitations and, thereby, to the ex- 
citon system. Both the selected vibrational mode and the 
exciton system will be described quantum mechanically, 
replacing the time dependence of the Hamiltonian (|l3) 
by a time-independent Hamiltonian that includes the se- 
lected vibrational mode and the exciton system. We as- 
sume that the electronic excitations of the B850 BChls 
form a linear system of coupled two-level systems, which 
interact at each site with dispersionless Einstein phonons 
of energy hiOQ. To simplify our notation, we will assume 
units with h = I. The corresponding model Hamilto- 
nian, the so-called Holstein Hamiltonian [p2i, reads 



H — Hex + Hph + Hint , 



iy^j 



Hph = ^ Wo b\ bi , 

i 

HM^guJoJ^BlB^ibl + b,) 



(60a) 
(60b) 

(60c) 
(60d) 



citation energy ef, Vij [c.f., Eq. (|lj)] encompasses the 
coupling between excitations of BChl^ and BChl^; Hph 

represents the selected vibrations of the BChls where b] 
and bi denote the familiar harmonic oscillator creation 
and annihilation operators. Hint describes the interac- 
tion between excitons and phonons, scaled by the dimen- 
sionless coupling constant g. 

The stationary states corresponding to the Hamilto- 
nian (pG) are excitons "dressed" with a phonon cloud and 
are referred to as polarons. Accordingly, the suggested 
description is called the polaron model. In what follows 
we shall restrict ourself only to the "symmetric case" , i.e.. 



-VS. 



Here Hex describes the electronic excitations of BChls; 
Bj and Bi are creation and annihilation operators ac- 
counting for the electronic excitation of BChl^ with ex- 



we assume e^ = eg for all sites, and we set Vij 
where V is the nearest neighbor interaction energy. This 
choice of the Hamiltonian is motivated by the results of 
our simulations, which suggest that the exciton dynamics 
is determined by the size- and time-scales of the thermal 
fluctuations of the excitation energies, e^, of individual 
BChls. Indeed, as shown in Fig. ||, the magnitude of the 
fluctuations of e^ is two orders of magnitude larger than 
the one corresponding to the nearest neighbor coupling 
energies, Vi,i±i, and at least three orders of magnitude 
larger than the one corresponding to the coupling ener- 
gies between non-neighboring BChls {Vij , with j ^ i ± 1) . 
Furthermore, according to Fig. 0, the period of oscilla- 
tion oiti(t) (i — \, . . . . 16) of a single BChl is about 20 fs, 
which corresponds to a high frequency intramolecular vi- 
bronic mode of energy o^o = 1,670 cm^^ = 207 mcV 
(see also Fig. p|). In the exciton Hamiltonian we set 
e^ = eo = 1.57 eV, and V = 350.9 cm'^ = 43.5 mcV, cor- 
responding to an exciton bandwidth of AV = 174 meV. 

The effect of static disorder can be incorporated into 
this model by deflning ti and Vij as random variables 
drawn from a suitably chosen distribution, which best 
characterizes the nature of the disorder. In this case, 
in calculating different physical observables, besides the 
usual thermal average, a second, conflguration average 
needs to be performed. 

The stationary states of the Holstein polaron Hamilto- 
nian (pOf) cannot be described analytically and one needs 
a suitable approximation. In general, well controlled per- 
turbative approximations are available only in the weak 
(5^1; for a more precise criteria see below) and strong 
(5^1) coupling limits. Besides conventional perturba- 
tion theory, the cumulant expansion method provides a 
convenient alternative for solving the polaron problem. 
The advantage of this method is that it provides reliable 
results even for arbitrary values of the coupling constant 
g [£^. By using our computer simulation results, we will 
estimate the value of g and show that our system falls 
into the weak exciton-phonon coupling regime. Then, we 
will use both perturbation theory and the cumulant ex- 
pansion method to investigate within the framework of 
the polaron model the effect of thermal fluctuations on 
the exciton bandwidth, coherence size, and absorption 
line-shape. 
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It should be noted that, because we are deahng with a 
single exciton coupled to a phonon bath, (i) the retarded 
Green's function will coincide with the real time Green's 
function Gfi(i) = G{t) = -i{TB{t)B^O)) (here T de- 
notes the time ordering operator), and (ii) the exciton- 
phonon coupling will not alter the phonon Green's func- 
tion. The latter can be written as iBSl 



Dit) = -i{TA{t)A{0)) 



'i[iNo + l)eM-i^o\t\) (64) 
-iVoexp(zcjo|t|)] , 



1.2 



1.4 



1.6 



1.8 



2.2 



where 



A{t) = b exp{~iujQt) + b^ exp{iujQt) 



is the phonon field operator, A^o = l/[exp(/3Li;o) — 1] is the 
Bose-Einstein distribution function, and (3 — l/fcsT. At 
room temperature (T w 300 K), we have f3ujo « 8.1, and 
consequently A^o ~ 3.1 x 10^"* <C 1. Thus, as we have 
already mentioned, the thermal population of the high 
frequency vibrational mode wq is negligibly small even at 
room temperature. 

In order to calculate the Green's function we sepa- 
rate the Hamiltonian (p^) into two contributions Ti = 
Ho + Hint , where Ho accounts for the first term in (|61 
Accordingly, we express 



8(eV) 

FIG. 9. Histogram of the fluctuations p{e) of individual 
BChl excitation energies e^ obtained from quantum chem- 
istry calculations. The vertical bars represent the weights 
of the Dirac-delta functions in the corresponding analytical 
result Q. 

A. Evaluation of coupling constant g 

The coupling constant g can be evaluated by estimat- 
ing the effect of thermal fluctuations on the electronic ex- 
citations of individual BChls. This corresponds to setting 
Vij = in (|6^). In fact, in calculating excitation energies 
ei{t) in our simulations we had not accounted for exci- 
tonic coupling. As a result, we can use ei{t), i = 1, . . . , M 
obtained from the simulations, to calculate the corre- 
sponding distribution function histogram [or density of 
state (DOS)] p{e) = dN{e)/de (see Fig. ||), as well as, any 

of the corresponding moments (e") = {l/N) J2j=i ^?i^j)- 

The coupling constant g can be determined by fitting 

to these data the prediction from Hamiltonian (p0|), for 

Vij = 0. The statistics for p{e) can be improved by av- ^^^^^ exp[-<i>(t)] = Go(t) exp 

eragmg over all sixteen B850 BChls. bmce we assume 

that all sites are equivalent, we can restrict ourself to the 

model Hamiltonian 



(65) 



G{t) = ^i{TB{t)B^O)) = Goit) cxp[-$(i)] , (66) 



where 



Go(t) = -ie(t)exp(-ieoO 



(67) 



is the Green's function corresponding to Ho- We then 
employ the cumulant expansion method [p3[ writing 



- E *™(*) 



m=0 

(68) 



H = eoB^B + ujob^b + guJoB''B{b + b'') 



(61) 



where 



where, for brevity, we dropped the irrelevant "i" site in- 
dex. 

The DOS p{e) can be obtained from the imaginary 
part of the Fourier transform of the retarded Green's 
function Gnit) = -ie{t){B{t)B'f{0)) [||, where (...) 
denotes the thermodynamic average over the vibrational 
(phonon) degrees of freedom in the exciton vacuum. Ac- 
cordingly, we employ the identity 



W„(i) = 



H? 



dti 



dU 



(2m)! 

•K{TB{t)nM{tl) . . .H^nt{t2v^)B\Q)) 



(69) 



P(^) 



1 



ImG_R(e) , 



where 



/oo 
dtGn{t)e''' ■ 
-OO 



(62) 



(63) 



and Hint{t) = guJoB^BA{t). The cumulants $m are ex- 
pressed in terms of Wm by identifying on both sides of 
Eq. (pq) all the terms which contain the same power of 
Hint- By using Eqs. (|6|,^) we have 

Wiit) = -{gu;o)^Goit) f dh f ' dt2D{h - ^2) (70) 

Jo Jo 
= Goit)[iept — So + S+ exp(— iwot) + -S*- e-xjp{iujot)] , 

where 

Cp = g'^ujo (71) 
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is the polaron binding energy, and 

So = g'(27Vo + 1) , ^± = g^No + 1/2 ± 1/2) . (72) 

From these results follows 

$i(i) = -G{7\i)W^i(t) = -iept + So - S+ exp{-iujot) 

—S- e'}q){iujot) . (73) 

By employing Wick's theorem |5^, one can calculate 
Wm, for arbitrary integer m > 1, with the result Wm = 
Go[Go Vi]"/to!. Consequently, we find $^(t) = for 
m > 1. The exact exciton Green's function is then given 
by 

G{t) = -ie{t) exp[-i(eo - ep)t - <i>o{t)] , (74a) 

^o{t) = So - S+ exp(-iwoi) - 5*- exp(iwoO • (74b) 



In order to calculate the DOS, we express Eq. ( |74b| ) in 
the equivalent form 



Mt) = e-a (2^0+1) exp{2gV^o(A^o + 1) (75) 

X cos,[iijJo{t + 1(3/2)]} , 



and employ the identity 

oc 

e =2^ h{z)i 



(76) 



where l(,{z) is the modified Bessel function. From 
Eqs. (|6^,|7^ [76|) one obtains finally the well known re- 
sult ||53,|54|| 



P(e) = X! Pi^i^-^d > 



(77a) 



with 



pt = exph.92(2iVo + 1) + £(/3c^o/2)] I,(2gV^o(A^o + 1)) , 

(77b) 

and 



= cjo^ + eo - ep 



(77c) 



In general, when the coupling Vij between excitons can- 
not be neglected, the higher order cumulants $„ do not 
vanish, and the Green's function G{t) cannot be calcu- 
lated exactly. Nevertheless, even in this case, the cumu- 
lant expansion method, as illustrated above in deriving 
the well known result (77), can be conveniently applied 
to calculate perturbatively, in a systematic way, both the 
Green's function and the corresponding absorption spec- 
trum (see below) to any desired degree of accuracy. 

The calculated DOS (normalized to unity) is an in- 
finite sum of Dirac-delta functions at energies e^, with 
weights Pi. As shown in Fig. ^, for a proper choice of 
the exciton-phonon coupling constant g, the "stick-DOS" 



(i.e., the weights pi) matches well the histogram p(e) ob- 
tained from our quantum chemistry calculations. The 
coupling constant g was determined by matching the first 
few moments of e, once calculated according to Eqs. (77a- 
77c), and then from the time series {ei(i)}i=i,...,i6- The 
definition of the moments 



/oo 
dep{e)e-=Y, 
-CX) « 



Pl^l 



(78) 



applied to (77a- 77c), results in analytical expressions 
that can be compared with the numerical values deter- 
mined from the time series. One obtains 

(e) = eo « 1.57 eV , (79a) 

(e2) = (e)2 + g^{2No + 1)cjo « 2.32 eV^ , (79b) 

(e^) = (e)3 + g^^l [l + ?,g\2No + 1) (79c) 



-K3(2iVo + 1) 



uo 



3.56 eV 



The coupling constant can be obtained from Eq. ([79q) 



9 = 



(62) 



(2iVo + l)uo'o 



0.65 , 



(80) 



Equation (79c) now contains no adjustable parameters 
and it can be used to check the reliability of our approach. 
It is found that the difference between the two sides of 
this equation is less than 2%. 

The value of 5, given by (pQ), does not determine by 
itself the value of the ratio between Hint (~ Cp) a-nd the 
problematic hopping term (~ AV, i.e., the energy band- 
width of the excitons) in the Holstein Hamiltonian. The 
actual dimensionless coupling strength parameter for the 
polaron model is k = ep/AV = g'^cuo/AV w 0.5. This 
value corresponds to a weak coupling regime of the po- 
laron model i3|,|5l. 



B. Polaron bandw^idth 

We return now to the full Holstein Hamiltonian (|6C 
with Vij 7^ 0. In the weak coupling limit it is convenient 
to rewrite the polaron Hamiltonian ( |60| ) in "momentum 
space" as 



H — Ho+ Hint , 

Ho = Y, ^kBlBk + wo ^ bib, 



Hint — 



? ' 

fc q 

k.q 



(81a) 
(81b) 

(81c) 



where Bk = {l/VM)J2jBjexp{ikj), etc. For Vij = 
—VSj,i±i {V > 0), the exciton dispersion (measured from 
the site excitation energy cq) is e^ = — 2ycos(fc), with 
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k = km — iTimjM ^ m = 0, 1,...,M— 1, and q = qp — 
27rp/M, p = 0,l,...,Af- 1. 

The polaron (phonon renormalized exciton) spectrum 
Ek can be calculated by using second order perturbation 
theory (all odd order terms vanish because the phonon 
creation and annihilation operators must appear in pairs) 



Ek = £k- 



9 ? 

M 



1 



e/c 



Efe+g — t^O 



(82) 



The renormalized exciton wave functions to leading order 
in g are given by 



\k;0) = \k;0) 






1 



e/c — e/c+q — ^0 



\k + q;lq) 



(83) 



Here |fc;0) = (1/vM) X^m ^''^P(*^"^)l"^) denotes the ex- 
citon wave function (in site \m) representation) in the ab- 
sence of the phonons, while \k + q; Ig) represents a state 
of non-interacting exciton (with momentum k + q) and 
phonon (with momentum q). The applicability of per- 
turbation theory is subject to the condition gujo/yM <^ 
|efc — Sk+q — Wo I, which in our case is marginally fulfilled 
since guJo/y/M — 33 meV< minjefc— efc+g— wol = 35 meV. 
According to Eq. ( p2|) the renormalized exciton (or po- 
laron) bandwidth is 



AEp = nia.x[Ek] ~ min[£'fe] w 66meV 



(84) 



which represents a 60% reduction with respect to the un- 
perturbed bandwidth AV. Thus the effect of the high fre- 
quency phonons is to reduce the exciton energy band, a 
situation commonly encountered in other weak coupling 
polaron models in which the phonons tend to enhance the 
quasiparticle mass, which is translated into a reduction 
of the width of the corresponding energy band |^, Q . 
However, it is known that static disorder (or coupling 
to low frequency phonons) has precisely the opposite ef- 
fect of increasing the exciton bandwidth in BChl aggre- 
gates. The fact that in our MD simulations we observe 
an increase instead of a decrease of the renormalized ex- 
citon bandwidth is most likely due to the fact that in 
our simplified polaron model the effect of static disorder 
is entirely neglected, while in the MD simulations this is 
implicitly included. 



C. Polaron coherence length 

Next, we determine the effect of the phonons on the 
coherence length Lp of the excitons in the BChl ring, by 
using our polaron model. At zero temperature, the exci- 
tons in the coupled BChl aggregate are completely delo- 
calized, and Lp should coincide with the system size M. 
In general, finite temperature and any kind of disorder 
will reduce the coherence length of the exciton. There is 



no universally accepted definition of Lp. We employ here 
an expression based on the concept of inverse participa- 
tion ratio, familiar from quantum localization p3|,|55[| 



Lp = 



E 



Pij\ 



^^E 



Pyl 



(85) 



where the reduced exciton density matrix is given by 



Py 



= J2cmCkij)expi-PEk) 



(86) 



In the absence of the exciton-phonon coupling, at room 
temperature, by setting in the above equations Ck{j) = 
ex-p{ikj)/y/M and Ek = e^ = —2Vcos{k), one obtains 
LpQ ~ 6.4, which represents a dramatic decrease with 
respect to the corresponding zero temperature value of 
16. For weak exciton-phonon coupling, the effect of the 
phonons on Lp can be taken into account via perturba- 
tion theory. By emp loying Eqs. (p3,p3[) in the density 
matrix (p6[), Eq. (|85| ) yields Lp w 5.4. As expected, the 
phonons, which act as scatterers for the excitons, reduce 
the coherence size of the latter. 

As it can be inferred from Fig. HG, the localization of 
the exciton is due primarily to thermal averaging and to 
a lesser extent to dynamic disorder. This conclusion is in 
agreement with previous studies |14[|. 



D. Absorption spectrum for single Einstein phonon 

Finally, we discuss the effect of the phonons on the 
absorption spectrum of the excitons. The fundamental 
absorption spectrum I{lu) (line shape function) is defined 
as [cf. Eq. (p9|)] 



I{lu) = 



1 

2^ 



dtf{t)e 



(87) 



where the generating function /(i), up to irrelevant fac- 
tors, is given by l54] 



/(i)..^|dfcpe-'=*((fc|C/(t)|fc)) 



(88) 



Here \dk\ is the magnitude of the transition dipole mo- 
ment connecting the ground electronic state and the 
|fc) exciton state; ((A:|...|fc)) denotes the thermal av- 
erage over the phonons of the excitonic matrix element 
{k\ . . .\k), and U{t) is an evolution operator given by 



U{t) = exp{iHot) exp[-i{Ho + Hint)t] 



(89) 
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FIG. 10. Top: Temperature dependence of the exciton co- 
herence length in the absence of disorder. Bottom: Polaron 
coherence length as a function of the exciton-phonon coupling 
constant g at room temperature. 



Hmtit) = exp(iiJof)iJ„f exp(-i_ffot) 



(90a) 



expressed in the interaction representation. The time 
dependent exciton density operator is given by 



Pq 



(i) = ^g-.(..-..+,)tS^S^ 



while the phonon field operator reads 



Mt) 



= bgc-'"^"' + blge"^°' 



(90b) 



(90c) 



The generating function (^) can be calculated by using 
the cumulant expansion method discussed above. Indeed, 
we can write in analogy to Eqs. (pl,p9h 



((fc|C/(i)|fc))=exph$fe(i)]=exp 



- E *™(^) 



rn—1 



(91) 



m=0 



where 



(-z)2 



(2m) 



dti 



dt 



2m 



(92) 



x{{k\TH,„tih) . . . H,„t{t2m)\k)) 



As for the Green's function (p^), Wm{t) contains only 
even numbers of time ordered Hint{t) oc A{t) operators, 
because the average over odd numbers of phonon field 
operators A(t) vanishes. Similarly to Eq. (vw, we have 



X{k\pg,{tl)pg^{t2)\k) 



2j2iTAgAtl)Ag,{t2)) 

91,92 



(5^0) 



\2 rt r-t 

7- / dh dt2y2'D{ti-t2) 

^ Jo Jo 



2M 

X{k\pg(tl)p^g{t2)\k) . 

A straightforward calculation yields 

Fkih -t2) = ^ J2(k\pg{h)p^g{t2)\k) 



M 

9 

-E 

9 



(93) 



(94) 



exp[i(efc - ek-g){ti - 12)] 



Inserting Eq. (M) into (|9^), and using the phonon 
Green's function (|64D, the cumulant ^i{t) can be written 



$i(t)=-VKi(i) = 



'f'LO^ i dT{t — T)i 

Jo 



D{T)Fk{T) . (95) 



If one neglects the coupling between the individual exci- 
tations (case corresponding to light absorption by indi- 
vidual BGhls), by setting T^ = 0, Eq. ( p5| ) yields the same 
expression as ([73|). Note that in this case, since tk = eo 
is k independent, we have F{t) — 1, and similarly to the 
calculation of the Green's function (|74|), it can be shown 
that W^{t) = [Wi{t)Y^/m\, for m > 2. As a result, 
all the corresponding higher order cumulants $m vanish, 
i.e., $(i) = $i(t), and the exact expression of the line- 
shape function /n(tj) for individual BChls assumes the 
familiar form i53,Ml 



/o(a;)ocexp(-5o) / dte*('^-^°+^'')* 



X exp(— 5+6 



— iu(]t 



-S- 



^iujQt 



E P£^(w - UJl) 



(96a) 



(96b) 



where pi and oji are given by Eqs. ( 77b ), and (77c), 
respectively. In other words, the absorption line-shape 
function of an exciton in the phonon field is given by the 



17 



imaginary part of the exciton Green's function, which is 
proportional to the density of states 1 54 1 . 

In general, /(w) cannot be calculated exactly. How- 
ever, even for arbitrary values of the exciton-phonon cou- 
pling constant 5, the cumulant approximation <i>fc(i) ~ 
^i{t) can be used safely to evaluate the generating func- 
tion f{t) and the corresponding line shape function I{lu). 
We have 



dT{t-T)V{T)Fk{T) 



Mt) 



where, according to Eq. 

V{t) = g^uj^iDit) 



iVne^"°*l 



The generating function is 



m - E 



dfcPe- 



-icfct -*fc(t) 



(97) 



(98) 



(99) 



and the corresponding line shape function is given by 
Eq. (p7|). Clearly, the coupling of the 16 exciton levels 
to a single Einstein phonon wq leads to a stick absorp- 
tion spectrum, i.e., a series of Dirac delta functions with 
different weights. 



E. Absorption spectrum for distribution of phonons 

In order to calculate, in the framework of the polaron 
model, the broadening of the absorption spectrum, and 
to compare it with the corresponding experimental result 
[see Fig. |8|, one needs to include the coupling of the ex- 
citons to the quasi-continuous distribution of the rest of 
the phonons. Formally, this can be achieved by replacing 
V{t) in Eq. IM) with [cf. Eq. iMj] 



V{t) 



where 



Jo 



it), (100) 



iD^{t) = {N^ + 1) cxp{-iujt) + N^ exp{iujt) . (101) 
Here we have introduced the phonon spectral function 

^H=E5X'5(c^-'^a). (102) 

a 

Once the actual form of Jjuj) is k nown, I{uj) can be cal- 
culated using Eqs. (^,p9|,p7p00|) . 

Apparently, the determination of J(w) requires the 
seemingly unattainable knowledge of the energies uja of 
all phonons, together with their corresponding coupling 
constants ga. However, the same problem posed itself 
in the framework of the spin-boson model description of 
the coupling between protein motion and electron trans- 
fer processes |M, and could be solved then through a 



spectral function evaluated from the energy gap fluctua- 
tions 6e(t). Likewise, here we can determine J{uj) from 
the fluctuations of the excitation energies 5e(i) calculated 
for individual BChls. The latter quantity can be obtained 
from the combined MD/quantum chemistry simulation 
carried out in this study. 

Indeed, the Hamiltonian for an individual BChl inter- 
acting with a phonon bath can be written 



H = Ho + Hrnt = (eo + 5€)B^B 



(103) 



where we defined the phonon induced energy-gap- 
fluctuation operator 



5e{t) = '^gaUJaAait) . 



(104) 



The autocorrelation function of the energy gap 5e{t) 
is the real part of the autocorrelation function of the 
energy-gap-fluctuation operator Si{t), i.e., 



C{t) = (Se{t)Se{0)) = Re[{Si{t)Si{0))] 



(105) 



Re 



9>UD^{t) 



Repit)] , 



wher e w e have use d Eqs. (104) and (100). Inserting 
Eq. (101) into (105), one obtains 



/■oc 

C{t) = duj J {uj)coth{f3uj/ 2) COS cut . 

Jo 



(106) 



J{lu) can be obtained through the inverse cosine trans- 
form of C{t), i.e., 



J(u;) 



tanh(/3w/2) / dtC{t) coswi 



(107) 



The autocorrelation function C(t) can be evaluated nu- 
merically from the flnite time series 5tj{ti). Here j = 
1,...,16 is the index labeling individual BChls, and 
ti — (i— 1) X 2 fs, z = 1, . . . , 400, denotes the time at which 
the energy gap was determined. As mentioned in the 
previous sections, each time series consisted of TV = 400 
time steps of 2 fs. For best sampling, one averages over 
all M = 16 BChls resulting in the time series 



M 



1 



N ■ 



N-i 

E 
fc=i 



5ej{ti + tk)6ej{tk) 



(108) 



whic h is shown in Fig. [ll|a. It is safe to assume that 
Eq. ( |l08| ) is reliable for t < 400 fs. We note that the au- 
tocorrelation function in Fig. [ll|a, resembles the transi- 
tion frequency autocorrelation function for Nile Blue [^ , 
obtained through a model that combines in both the in- 
tramolecular vibrations and solvent dependent contribu- 
tions. 

The corresponding phonon spectral function J(a;) can 
be computed by employing Eq. (107). The result is shown 
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FIG. 11. (a) Autocorrelation function C{t) of the energy 
gap fluct uati ons Se{t) for individual BChls, calculated us- 
ing Eq. (lOS). The inset shows the short time behavior of 



C{t). ( h) P honon spectral function J{uj) obtained according 
to Eq. (p^). 



in Fig. [Tllb. The shape of J(a;) resembles the power spec- 
trum of the BChl excitation energies (see Fig. @). In 
particular, the prominent peak around loq ~ 0.2 eV in 
J(w) indicates a strong coupling of the system to an in- 
tramolecular C=0 vibronic mode. Note, however, that 
J {to) has significant contributions over the entire range 
of phonon energies < w ^ 0.22 eV and, therefore, we 
expect that all, i.e., low, intermediate and high frequency 
phonons will contribute to the broadening of the line 
shape function I{oj). 

By assuming that the coupling of BChl molecules to 
phonons (intramolecular vibrations, vibrations of the 
protein matrix and solvent molecules) is independent 
of their excitonic coupling, we may conclude that J{uj) 
given by Eq. (107), describes equally well both individ- 
ual BChls and excitonic aggregates of BChls. The only 
difference between the corresponding absorption spectra 
comes from the "excitonic" factor Fk(t) in Eq. (|97|). 

In order to calculate the line shape function I{uj), as 
given by Eqs. (^,^,^, from the energy gap autocor- 



relation function C{t) one proceeds as follows. First, we 
determine 

V{t)=Vi{t)-iV2{t) (109) 

dLjjJ{uj)[coi\i[j3Lo/2) cosLot — isvuut] , 



from which, by taking into account Eqs. (106) and (107) 
one obtains 



and 



I?i(i)=C(i), 



T^2{t) — / du!j{uj)smLjt 



(110a) 



(110b) 



Next, we calculate the cumulant ^k{t)', according to 

Eq.m 



$fe(i)-$'fc(i)-«*"fcW, 



(lUa) 



with 



$'fc(t) = / dr{t - r)[C{T)Re{Fk{t)} - I?2(r)Im{^fc(i)}] 



(111b) 



and 



$"fe(i) = / dT{t - T)[D2{T)Re{Fk{T)}+C{T)Im{Fk{T)}] 
Jo 

(111c) 

Finally, the line shape function is according to 
Eqs. (g,|9D 

/•OO 

/(cj)(x V|dfc|2 / diexp[-$'fc(i)]cos[(a;-efc)i + $"fe(t)] 
k -^0 

(112a) 

This result applies to the general case when the system 
has several optically active levels, characterized by differ- 
ent transition dipole moments dk- In both cases consid- 
ered by us, i.e., individual BChls and the fully symmetric 
B850 system (i.e., a circular aggregate of 16 BChls with 
excitonic coupling, having 8-fold symmetry, and transi- 
tion dipole moments oriented in the plane of the ring of 
BChls), there is i n fact only one optically active level and, 
therefore, in Eq. ( 112a ) the summation over k, as well as 
the transition dipole moment can be both dropped. For 
an individual BChl we take efc = eq ~ 1-6 eV, while for 
the B850 system the optically active doubly degenerate 
level is Cfe = e±i — eo — 2I/cos(7r/8) — 1.52 eV. In this 
case, the line shape function assumes the simpler form 

/•CJO 

I{uj) cx / dtexp[-$'fc(t)] cos[{uj - ek)t + $"fe(i)] • 

(112b) 
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FIG. 12. (a) Absorption spectra of individual BChls: I{ui) 
- time series analysis (dashed curve), Io{oj) - polaron model 
with <E>"o = (thin solid curve), iBChiis^) - polaron model 
(thick solid curve), (b) Comparison between the normalized 
DOS p{i^) and Ischiis^), obtained within the framework of 
the polaron model. 



In case of individual BChls, i.e., without excitonic cou- 
pling, F};{t) — 1, and according to Eqs. (112b, rTTh| ,lllc) 
the absorption spectrum reads 



lBChi{uj) oc / diexp[-$'o(i)] cos[(w - eo)t + $"o(<)] , 

(113a) 



where 



and 



$'o(t) - / drit - t)C{t) , 
Jo 



$"o(t)= / dr(t-r)I?2(r). 



(113b) 



(113c) 



For the moment let us neglec t the imaginary part of the 
cumulant $o(0 i^i Eq. ( 113a ), i.e., we set ^"o{t) w 0. It 
follows 



^BChii^) ~ h{^) « / dtexp[-$'o(i)] cos[(w - eo)i] , 

(114) 



and one can easily see that Iq{lo) coincides (up to an 
irrelevant normalization factor) with the cumulant ap- 
proximation of the line shape function /(w), Eq. ([5S|), 
derived and employed in our computer simulation study. 
Therefore, it comes to no surprise that in Fig. n3a the 
plots of the normalized [with max{/(cL;)} = 1] line shape 
functions Iq{uj) and /(w) completely overlap in the peak 
region. The difference between the curves away from the 
central peak is most likely due to poor statistics and nu- 
merical artifacts introduced by the FFT evaluation of the 
integrals for the time series analysis result (p9) . The peak 
of the absorption spectrum is located at eo — 1.6 eV and 
the corresponding FWHM is 125 meV. 

However, t he ac tual line shape function IbcmI'-^), as 
given by Eq. (113a), include a non- vanishing $"o(i). The 

As one 



12 i. 



corresponding result is also plotted in Fig. 
can see, the contribution of <I>"o(i) redshifts the peak of 
Inchii^) to 1.53 eV, and renders Ischii^) asymmetric 
with a somewhat larger FWHM of 154 meV. As illus- 
trated in Fig. p^, iBC hi(^) matches the corresponding 
DOS p{uj) [see Eq. (77a) and Fig. y. This is not surpris- 
ing since, as we have already mentioned, for a two level 
system which is linearly coupled to a phonon bath, both 
the DOS and the absorption spectrum are proportional 
to the irnaginary part of the corresponding Green's func- 
tion [^3| , |54[| . The line shape function (|59| ) of individual 
BChls derived within the theoretical framework of Sec. II 
does not account for the imaginary part of the cumulant 
$o(i) which suggests that the time development operator 
C/(t,io) [Eq. (P0|)] is not totally adequate for calculating 
the absorption coefficient. The situation seems to be sim- 
ilar to the commonly used stochastic models in modeling 
the absorption spectrum of excitonic systems coupled to 
a heat bath with a finite time scale [£8|. In this case as 
well, the autocorrelation function of the stochastic en- 
ergy gap fluctuations is assumed to be real (i.e., $"(i) is 
neglected), and as a result the model violates the fluctu- 
ation dissipation theorem, and yields no Stoke shift [ p8[ . 
This discrepancy is due to the fact that, unlike in the case 
of the polaron model, the system-heat bath interaction 
is not incorporated in a consistent manner. 

Within the framework of the polaron model, the ab- 
sorption spectrum /_b85o(^) of the excitonically cou- 
pled B850 BChls can be calculated numerically by using 
Eqs. ( |1 1 2b| , [ril| , p0^JT0^ ,|9^) . The needed input quanti- 
ties are: T = 300 K, C{t), and e^ = e±i = 1.52 eV; these 
are derived from our simulations and there are no free 
parameters in the above mentioned equations. The com- 
puted /s85o('-^) is shown in Fig. O, along with the corre- 
sponding spectrum Itsa{^) obtained through time series 
analysis, [c.f. Eqs. ( p6| , p7|)] , and the corresponding exper- 
imental result /exp(t^) ||42|]. The peak of /ssso is located 
at 1.46 eV (w 849 nm), and coincides with the position 
of the other two spectra. /_b85o(w) (FWHM=88 meV) is 
somewhat broader than Iexp{i^) (FWHM=58 meV), but 
is narrower than /tsa(w) (FWHM=138 meV). The effect 
of exchange narrowing due to the excitonic coupling be- 
tween B850 BChls is manifest: FWHM is reduced from 
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FIG. 13. Normalized absorption spectra of the B850 BChls 
with excitonic coupling: Ibs5o{i^) ~ polaron model (solid 
curve), IrsAii^) - time series analysis (dashed curve) and 
lexpi'jj) - experiment (circles). 



154 meV (Fig. |l3|) to 88 meV (Fig. ^), which represents 
a reduction in the width of the line shape function by a 
factor of 1.75. This exchange narrowing is accounted for 
by the factor Fk{t) in Eq. (|5); since \Fk(t)\ < 1, this 
factor reduces the value of the cumulant <&fc(i), which 
leads to a narrowing of the line shape function [recall 
that $fc = result in a Dirac-delta function type /(y)]. 
In contrast to the commonly used stochastic models pq] 
for describing the absorption spectrum of molecular ag- 
gregates, our polaron model docs not postulate the form 
of the autocorrelation function C{t), but rather uses this 
function as an input, derived from computer simulations 
or from experiment. In this regard, our polaron model 
provides a more realistic approach for evaluating optical 
properties of molecular aggregates in general, and the 
B850 excitons in particular. Indeed, a generic autocor- 
relation function C{t) = A^exp(— At), assumed by most 
stochastic models, would represent a gross oversimplifica- 
tion of the complicated structure oiC{t) as inferred from 
our computer simulations (see Fig. [ll]a). It is clear that 
this autocorrelation function cannot be modeled by either 
a single or a sum of several exponentially decaying func- 
tions, in spite of the fact that by a proper choice of the 
values of the mean square fluctuations of the energy gap, 
A^, and of the inverse relaxation time. A, the stochastic 
model can yield an almost perfect fit to the experimental 
spectrum. For example, by assuming that in our case the 
broadening of the line shape function is due to coupling to 
phonons with clearly separated high and low frequency 
components, it can be shown that the real part of the 
cumulant is approximately $'(t) « A^i^/2 -I- Tht, where 
£ (h) refers to the low (high) frequency component and 
Tfi — Af/Xh. (The imaginary part of the cumulant $"(i) 
is assumed to affect only the shift of the absorption spec- 



trum, but not its broadening.) The first (second) term in 
^'(i) represents the short (long) time approximation of 
the cumulant brought about by the low (high) frequency 
phonons. The corresponding line shape function has a 
Voigt profile |5q], i.e., is a convolution of a Gaussian and 
Lorentzian due to low and high frequency phonons, re- 
spectively. For As = 25 meV and F^ = 7 meV, one can 
obtain an almost perfect fit to the experimental B850 
absorption spectrum. However, in spite of this apparent 
agreement with experiment, neither these numerical val- 
ues nor the clear separation in the phonon spectral den- 
sity is justified in view of our computer simulation results, 
and as such this kind of analysis should be avoided. In 
fact, as it can be inferred from the plot of J(a;) shown in 
Fig. [ll|b, all phonon frequencies up to 0.24 eV contribute 
to the cumulant <I>(i) and, therefore, to the absorption 
spectrum I {to). Thus, unlike stochastic models which 
employ fitting parameters to account for the broaden- 
ing of the experimental absorption spectra, our polaron 
model incorporates dynamic disorder directly through 
the phonon spectral function J(w) obtained from com- 
bined MD/quantum chemistry calculations, i.e., without 
any ad hoc assumptions and fitting parameters. 

Further tests of our polaron model based on the calcu- 
lation of other optical properties of the system, such as 
the circular dichroism (CD) spectrum, are in preparation. 



VI. DISCUSSION 

This paper presents a novel approach to the study 
of excitons in light-harvesting complexes that combines 
molecular dynamics and quantum chemistry calculations 
with time-dependent effective Hamiltonian as well as po- 
laron model type of analyses. 

The molecular dynamics approach allowed us to ob- 
serve, at the atomic level, the dynamics of BChl motions 
and gain insight into the extent and timescales of geo- 
metrical deformations of pigment and protein residues at 
physiological temperatures. A comparison of the average 
distances and angles between the B850 BChls to those 
found in the crystal structure reveals an increased dimer- 
ization within the B850 ring, as compared to the crys- 
tal structure. We observed diffusion of a water molecule 
into the B800 BChl binding site for seven out of eight 
B800 BChls. The average location of this water molecule 
agrees well with the crystal structure location. Future 
quantum chemistry calculations will determine whether 
this molecule plays a functional role. In this respect it is 
interesting to note that the recently published structure 
of the cyanobacterial photosystem I at 2.5 Aresolution re- 
vealed also ligations of Chla Mg ions involving side groups 
other than His, and, in particular, water p9| ]. 

The time series analysis theory, based on the combined 
MD/QC calculations resulted in an absorption spectrum 
that is about a factor of two wider than the experimen- 
tal spectrum. One of the reasons for this discrepancy 
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might be improper treatment, and subsequently an over- 
estimate, of the contribution of high-frequency modes 
within the framework of our combined MD/QC calcu- 
lations. We also note that the time series analysis ne- 
glects some important quantum effects which may cause 
the discrepancy between the computed and the experi- 
mental results. These quantum effects are accounted for 
by the polaron model. The amplitude of fluctuations of 
the off-diagonal matrix elements, i.e., couplings between 
BChls, was found to be at least two orders of magnitude 
smaller than the corresponding fluctuation amplitude of 
the diagonal matrix elements. We believe that in spite 
of the possible overestimate of the fluctuation associated 
with the high frequency modes, we can safely conclude 
that the disorder in the B850 system is diagonal rather 
than off-diagonal in nature. 

The time-dependent effective Hamiltonian description 
revealed an exciton spectrum red-shifted as compared to 
the spectrum of individual BChls. This red-shift is well 
known in J-aggregates, and is attributed to the transfer- 
ring of the dipole strength into low-energy exciton states. 

The observation that the fluctuations of the BChl 
site energies could be attributed largely to high fre- 
quency intramolecular vibrational modes, with energy 
LUo r^ 1,600 — 1,700 cm^^, prompted us to model the 
effect of dynamic disorder on the B850 excitons by em- 
ploying a polaron Hamiltonian, which describes both ex- 
citons and the coupled single phonon mode by quantum 
mechanics. The strength of the exciton-phonon coupling 
g is related to the RMSD of the fluctuating site energies, 
and was found to be weak. By employing standard per- 
turbation theory we investigated the effect of dynamic 
disorder (i.e., exciton-phonon coupling) on the energy 
dispersion and localization of the excitons. We found 
that in contrast to static disorder, which leads to a broad- 
ening of the exciton bandwidth, dynamic disorder due to 
high frequency vibrational modes leads to a reduction of 
the width of the exciton energy band. Also, our polaron 
calculations showed that dynamic disorder reduces only 
slightly the exciton delocalization length (from around 
6 BChls to 5 BChls at room temperature), confirming 
previous results according to which the main mechanism 
responsible for exciton localization in LH-II rings is ther- 
mal averaging W4\. 

By employing the cumulant expansion method, we cal- 
culated, in the framework of the polaron model, the ab- 
sorption spectrum of the B850 BChl, with and without 
exciton coupling, for a single phonon mode as well as for 
a distribution of phonons. The absorption spectrum of 
the excitonic system coupled to a single high frequency 
phonon ujq (intramolecular vibronic mode) is given by a 
series of Dirac-delta functions (stick spectrum) , with dif- 
ferent weights. In order to obtain a realistic, broadened 
absorption spectrum, which can be compared with the 
one measured experimentally we had to included the ef- 
fect of the entire phonon distribution through the phonon 
spectral density J{ui). We have shown that J(w), which 
accounts for both coupling strengths and frequencies of 



the individual phonon components, can be determined 
from the autocorrelation function C{t) of the energy gap 
fluctuations of individual BChls, readily available from 
our combined MD/quantum chemistry simulation. We 
found that in our model only three inputs, namely au- 
tocorrelation function C(i), temperature T and optically 
active exciton level e^, are needed to calculate the absorp- 
tion spectrum. Considering the fact that these inputs 
come directly from our MD/QC calculations, and that 
we have no free parameters, the agreement between the 
spectrum obtained from the polaron model calculations 
and the experimental spectrum is remarkable. 
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